import numpy as np


def linear_regression(X, y):
    X_b = np.c_[np.ones((X.shape[0], 1)), X]
    X_b_T_X_b = np.linalg.inv(np.dot(X_b.T, X_b))
    theta_best = np.dot(np.dot(X_b_T_X_b, X_b.T), y)
    return theta_best